Water‐limited environments affect the association between functional diversity and forest productivity

Abstract The link between biodiversity and ecosystem function can depend on environmental conditions. This contingency can impede our ability to predict how biodiversity‐ecosystem function (BEF) relationships will respond to future environmental change, causing a clear need to explore the processes underlying shifts in BEF relationships across large spatial scales and broad environmental gradients. We compiled a dataset on five functional traits (maximum height, wood density, specific leaf area [SLA], seed size, and xylem vulnerability to embolism [P50]), covering 78%–90% of the tree species in the National Forest Inventory from Italy, to test (i) how a water limitation gradient shapes the functional composition and diversity of forests, (ii) how functional composition and diversity of trees relate to forest annual increment via mass ratio and complementarity effects, and (iii) how the relationship between functional diversity and annual increment varies between Mediterranean and temperate climate regions. Functional composition varied with water limitation; tree communities tended to have more conservative traits in sites with higher levels of water limitation. The response of functional diversity differed among traits and climatic regions but among temperate forest plots, we found a consistent increase of functional diversity with water limitation. Tree diversity was positively associated with annual increment of Italian forests through a combination of mass ratio and niche complementarity effects, but the relative importance of these effects depended on the trait and range of climate considered. Specifically, niche complementarity effects were more strongly associated with annual increment in the Mediterranean compared to temperate forests. Synthesis: Overall, our results suggest that biodiversity mediates forest annual increment under water‐limited conditions by promoting beneficial interactions between species and complementarity in resource use. Our work highlights the importance of conserving functional diversity for future forest management to maintain forest annual increment under the expected increase in intensity and frequency of drought.


| INTRODUC TI ON
A compelling amount of empirical evidence has demonstrated that biodiversity can affect ecosystem functioning, with many studies supporting positive biodiversity-ecosystem function (BEF) relationships (Cardinale et al., 2011;Chapin et al., 1997;Loreau et al., 2001;Ratcliffe et al., 2016;Wright et al., 2021). The importance of BEF effects can be context-dependent, changes in environmental conditions may alter the strength and shape of these relationships (Pretzsch et al., 2013;Ratcliffe et al., 2016;Wright et al., 2021). This impedes our ability to predict how BEF relationships will respond to rapid global change. As such, there is a clear need to explore the processes underlying shifts in BEF relationships across large spatial scales and broad environmental gradients (Baert et al., 2018;Fei et al., 2018;Gonzalez et al., 2020;Grossiord, 2020;Hisano et al., 2018;Paquette et al., 2018;Ratcliffe et al., 2017).
Characterizing communities in terms of their functional composition and diversity can help reveal the factors shaping patterns of biodiversity and community structure, as well as the influence of biological communities on ecosystem functioning (Augusto & Boča, 2022;Bonilla-Valencia et al., 2022;Dıáz & Cabido, 2001;Laureto et al., 2015;Petchey et al., 2004). Numerous studies have shown that considering the functional traits of organisms can improve our understanding of how biodiversity affects ecosystemscale processes by providing a more physiological basis for the ecological function of species in communities (Cadotte et al., 2011;Dıáz & Cabido, 2001;Song et al., 2014;Yan et al., 2023). For example, Ayma-Romay et al. (2021) showed that the inclusion of traits, which capture key variation in plant life-history strategies, representing a trade-off between conservative and acquisitive resource strategies, can help explore the mechanisms that underlie BEF relationships.
Two main mechanisms are generally considered to underlie BEF relationships: niche complementarity (or resource partitioning) and mass ratio (or dominance) effects (Loreau et al., 2001). Niche complementarity can be reflected by the functional diversity of a community; a greater variety of trait values in a given community can lead to higher rates of ecosystem function due to more efficient exploitation of resources (Sonkoly et al., 2019). Mass ratio effects refer to the composition of a community; the trait values of the most abundant species in a community are expected to have the largest effect on the relationship between diversity and ecosystem function (Ali, 2015;Loreau et al., 2001;Sonkoly et al., 2019). Previous studies have produced contrasting results regarding the relative importance of these two mechanisms to explain BEF relationships (Ammer, 2019;Grossman et al., 2018;Wright et al., 2021). However some common patterns have emerged; specifically, environmental stress (e.g., water limitation) appears to be a major factor governing the relative importance of niche complementarity versus mass ratio effects. For example, Richardson et al. (2012) showed that the relative importance of complementarity and mass ratio effects varied inversely with a latitudinal gradient in grasslands, with a greater relative importance of complementarity effects in more water-limited environments.
The stress gradient hypothesis (SGH; Bertness & Callaway, 1994) provides a framework for predicting how environmental conditions may mediate BEF relationships. Specifically, the SGH predicts that the frequency of facilitative and competitive interactions will vary inversely across stress gradients (Bertness & Callaway, 1994;Maestre et al., 2009). In relatively benign conditions (e.g., mesic forests), competitively dominant species will drive rates of ecosystem function via mass ratio effects. Under higher abiotic stress (e.g., xeric forests), in contrast, facilitation among species will lead to more efficient use of available resources at the community level and outweigh the negative competitive impacts of neighbors. As a result, the role of competitive interactions becomes relatively weak compared to niche complementarity effects via resource partitioning (Richardson et al., 2012;Schmitt et al., 2020).
To better understand how the processes influencing BEF relationships change across broad environmental gradients, we used a functional trait-based approach over an extensive spatial scale (i.e., Italy) with a strong climate gradient. We used a National Forest Inventory dataset to examine the mechanisms by which tree diversity influences annual increment and how the relative importance of complementarity and mass ratio effects changes along a broad gradient of water limitation. We addressed the following questions:

| How does the functional composition and diversity of forests vary with respect to abiotic gradients throughout Italy?
We expect community-weighted mean (CWM) trait values to vary with respect to abiotic conditions in a way that reflects selection across environmental gradients. Specifically, we expect more resource-conservative traits under increased water limitation ( Figure 1, Table 1). Additionally, in more water-limited environments, we expect lower functional diversity as a result of selection for a narrow range of more conservative growth strategies (Wieczynski et al., 2019; Figure 1). On the contrary, if competitive exclusion is the main driver behind functional diversity under favorable environmental conditions, we expect to see lower functional diversity under favorable environments in comparison with water-limited environments ( Figure 1; Levine & HilleRisLambers, 2010;Mayfield & Levine, 2010;Mensah et al., 2018).

T A X O N O M Y C L A S S I F I C A T I O N
Biodiversity ecology, Community ecology, Ecosystem services studies, Functional ecology

| How does functional composition and diversity of forests relate to annual increment?
We expect a positive relationship between tree diversity and annual increment, driven by a combination of niche complementarity and mass ratio effects (Figure 1; Ali, 2015;Loreau et al., 2001;Sonkoly et al., 2019). More specifically, we expect niche complementarity (indicated by functional dispersion) to be positively associated with site annual increment, whereas we expect the direction of mass ratio effects (indicated by CWM values) to vary with different traits (Conti F I G U R E 1 Illustration of the main research questions tested in this work. RQ1: Hypothesized response of functional composition and functional diversity against a water limitation gradient. Green triangles display increasing values, red triangles display decreasing values under increased water limitation. RQ2: Theoretical model of climate (VPD) and community functional properties (CWM and FDis) on annual volume increment (Cai). Solid arrows indicate predicted causal relationships among variables and lowercase letters are path estimates. RQ3: Stress gradient hypothesis provides a framework for understanding the mechanisms by which biodiversity influences annual volume increment, predicting that the frequency of facilitative and competitive interactions will vary inversely across abiotic stress gradients (i.e., water limitation gradient).

RQ2 RQ3
Water limitation gradient (VPD)  & Díaz, 2013;Finegan et al., 2015;Mensah et al., 2018). In general, we hypothesize that the dominance (i.e., mass ratio effects) of species with more acquisitive functional traits (e.g., higher values of height, SLA, and P 50 , lower values of wood density, and seed mass) to be positively associated with site annual increment.

| Is functional diversity more strongly related to annual increment in the Mediterranean or temperate climate region?
In general, forests in the Mediterranean climate region are ex-

| Study area and forest inventory data
The study area extends throughout Italy (35°29′-47°04′ N, 6°37′-18°31′ E, Figure S1), covering a highly variable climate gradient that ranges from Mediterranean to temperate climatic regions. More than one-third of the country's 30 million hectares of land area is covered in forests and other woodlands Gasparini et al., 2022). Oak-, beech-, and chestnut-dominated forests are the most common forest types, each representing over 10% package (Hijmans, 2020) in R v. 4.1.3 ( Figure S1; R Core Team, 2021).
We focus on vapor pressure deficit (VPD, kPa) to represent a gradient of water limitation as it is an integrative measure of water stress that reflects the effects of both precipitation and temperature.

| Functional trait data
To characterize the functional diversity of trees in INFC plots, we compiled publicly available data from the TRY (Kattge et al., 2020) and WOODIV (Monnet et al., 2021) databases. Note that a caveat of publicly available data is that we lack trait data for some rare species, however, the coverage of the trait data was high (see Appendix S1).
We focused on the following traits: maximum tree height (m), seed mass (mg), wood density (g cm −3 ), specific leaf area (SLA; mm 2 g −1 ), and xylem vulnerability to embolism measured as the xylem pressure inducing 50% loss of hydraulic conductivity due to embolism, that is, P 50 (MPa) ( Table S1). We chose these traits because they capture For each trait (including the multivariate trait axis, PC1), we calculated the plot-level community-weighted mean value (i.e., functional composition, CWM) and a measure of functional diversity (i.e., functional dispersion, FDis) using the R package "FD" (Laliberté & Legendre, 2010). We used species relative abundance to weigh both indices. FDis is generally independent of species richness (in our dataset, Pearson's r ranged from .04 to .15, depending on the trait; Figure S3).

| Statistical analyses
Separately for each trait (including PC1), we used Structural Equation Modeling (SEM) to test mediation hypotheses for relationships among community functional properties (composition and dispersion) and annual increment. Our conceptual a priori model (Figure 1) is based on previous research designed to disentangle the relationships among ecosystem function, functional diversity, and community-weighted mean traits (Chiang et al., 2016;Lohbeck et al., 2015). We applied structural equation modeling to the entire dataset, compiling a series of structural paths for each trait (i.e., seed mass, maximum height, SLA, wood density, P 50 , and PC1). Prior to analysis, all variables were log-transformed (i.e., log(x + (1 − min(x)))), centered, and scaled to standard deviation to account for negative values, improve the symmetry of the distributions, and facilitate model fitting. A Wishart likelihood approach was used for the maximum likelihood (ML) estimation, and a full information maximum likelihood (FIML) method was used for missing data (Wothke, 1998 Mediterranean. This multigroup analysis provides a direct test of measurement invariance between climatic groups, thus ensuring that the observed differences in structural relationships across conditions are unaffected by neither measurement errors nor measurement differences (see Appendix S1). SEMs were conducted using the "lavaan" package (Rosseel, 2012) implemented in the R environment v 4.1.2 (R Core Team, 2021).

| RE SULTS
Across Italy, the water limitation gradient (i.e., VPD) ranged from 0.14 to 0.90 kPa, while in the temperate climatic region it ranged from 0.14 to 0.78 kPa (mean 0.47 ± SD 0.12) and in the Mediterranean climate region it ranged from 0.36 to 0.90 kPa (mean 0.61 ± SD 0.09).
Functional traits varied considerably across species (Table S1); PC1 explained 36.1% of the total trait variation and primarily represented a trade-off between conservative and acquisitive resource strategies (i.e., negative correlation with SLA and height; Figure S2, Table 1). Forest annual increment varied substantially across Italy, ranging from 0.005 to 13.672 m 3 ha −1 year −1 .

| Functional composition and dispersion along a climate gradient
The fit indices of specified SEM and MG-SEM models, the former based on pooled data ( Figure 2) and the latter on grouped data ( Figures S4-S9), are well within the acceptable limits ( Figure 2, Figures S4-S9). Patterns of functional composition and dispersion across the VPD gradient depended on spatial scale and varied between the climate regions (Figures 3 and 4, Tables S3-S14). At the national scale, CWM was positively associated with VPD for seed mass, SLA, wood density, and PC1 but negatively associated with VPD for maximum height. These trends were consistent when considering only the temperate forest plots and, in addition, CWM of P 50 increased with VPD. Among plots in the Mediterranean climate region, VPD was negatively associated with CWM for seed mass, maximum height, SLA, and P 50 and positively associated with CWM for wood density (Figure 3, Tables S3-S14). At the national scale, FDis was positively associated with VPD for seed mass and P 50 but negatively associated with VPD for SLA. These trends were consistent for the temperate forest plots and, in addition, FDis of maximum height, wood density, and PC1 increased with VPD. Among plots in the Mediterranean climate region, VPD was negatively associated with FDis for wood density and positively associated with FDis for P 50 (Figure 4, Tables S3-S14).

| The link between functional composition, diversity, and annual increment
Paths linking annual increment to selected variables (i.e., VPD, CWM, and FDis) explained a total variance ranging from 7% to 10% based on pooled data and from 4% to 12% based on data grouped by climate region (Tables S3-S14). All parameter estimates and related fit indices are shown in Tables S3-S14. Our separate SEM models for each trait reveal that both CWM values and FDis were significantly related to annual increment for all traits, with the exceptions of maximum height, which was significantly related to annual increment only through CWM, and PC1, which was only significantly related to annual increment through

| Functional composition and diversity along a climate gradient
Although specific trait shifts along environmental gradients have commonly been observed (Costa-Saura et al., 2016, 2019Joswig et al., 2022;Pinho et al., 2021), we only found a consistent direction of functional composition shifts in the two climate regions for three traits (i.e., height, wood density, and PC1), highlighting the importance of spatial scale when assessing plant trait response to abiotic factors. In our study, opposite ends of the water limitation gradient likely face differing limiting factors. In the temperate climatic region, cold temperatures limit plant growth and development (Körner et al., 2016), and freezing stress can lead to mortality (Pittermann & Sperry, 2006). In the Mediterranean region, however, water limitation is a major limiting factor for plant growth and development (Gazol et al., 2018). The combined effect of cold stress and water limitation across the entire VPD gradient could help explain why we found a hump-shaped response for certain traits (i.e., SLA, Seed mass, P 50 ), with a shift toward coniferous species in cold stressed mountainous areas being the key factor behind more conservative trait values under low VPD in the temperate regions (Charrier et al., 2021;McCulloh et al., 2023).
Nevertheless, in the Mediterranean climatic region, we found that CWM traits tended to be more conservative in sites with higher levels of water limitation (higher VPD). As hypothesized, we found a trend of decreasing SLA, height, and P 50 , while stem den-  Larios and Venable (2018) showed that under water-limited conditions when the cost of construction is considered, there is no overall fitness increase with seed size.
As with trends of CWM traits, variation of FDis along the VPD gradient also depended on the trait considered, the spatial scale of analysis, and the climate region of the plots. Nevertheless, we found a consistent direction of FDis shifts among the temperate plots, FDis was positively associated with VPD (i.e., higher functional diversity in more water-limited sites) for seed mass, height, wood density, P 50, and PC1. At the same time, FDis of SLA was lower in sites with higher VPD.
FDis of SLA is expected to be highest when coniferous and angiosperm species co-occur due to large differences in leaf strategies (Maynard et al., 2022). This suggests that an increasing VPD is causing a shift toward the dominance of angiosperm species (Charrier et al., 2021;McCulloh et al., 2023). The FDis pattern of SLA along the VPD gradient confirms that the hump-shaped trends for CWM traits, such as SLA, seed mass, and P 50 , result from a shift between coniferous and angiosperm species due to cold stress at lower VPD values. Additionally, an increase in FDis among the other five traits indicates a broad range of functional strategies that is able to co-occur. The release of cold stress with increasing VPD in the temperate region leads to relatively benign environmental conditions, resulting in increased functional diversity.

| Links between functional composition, diversity, and annual increment
Overall, our results support a positive relationship between functional diversity and annual increment (Cardinale et al., 2011;Chapin et al., 1997;Loreau et al., 2001). More specifically, we found evidence for niche complementarity effects through a positive association between FDis and site annual increment, suggesting that resource partitioning within tree communities positively influences annual increment in Italian forests. In contrast, maximum tree height was only associated with annual increment through mass ratio effects (CWM value), which aligns with previous research in that mass ratio effect of maximum height is strongly related with forest annual increment (Chiang et al., 2016;Conti & Díaz, 2013;Finegan et al., 2015). On the contrary, the direction of mass ratio effects varied among the traits; sites dominated by species with more acquisitive traits (e.g., increased height, SLA, and P 50 ) showed higher annual increment, on average. Our results suggest that annual increment of Italian forests is positively influenced by a diversity of resource strategies that allow for resource partitioning, while at the same being positively associated with the dominance of tree species with more acquisitive resource strategies.

F I G U R E 3
Results of Structural equation modeling (SEM) on the pooled dataset (black) and multigroup structural equation modeling (MG-SEM) for the temperate (blue) and Mediterranean (red) bioclimatic domains. Lines represent the paths between community-weighted mean (CWM) of seed mass (SeedMass), tree height (Height), specific leaf area (SLA), wood density (WD), xylem vulnerability (Xylem) function traits, and all traits (All) and a vapor pressure deficit (VPD); solid lines represent significant (p < .05) paths, dashed lines not significant ones.  Drought is expected to increase in intensity and frequency in Italian forests (Spinoni et al., 2018), which could further alter the composition and annual increment of these forests. On one hand, we showed that those tree communities tended to be more conservative in sites with higher water limitation, while on the other hand, annual increment was positively influenced by the dominance of tree species with more acquisitive strategies, suggesting a decrease in forest annual increment in sites with higher water limitation due to a shift in functional composition. Under increased drought conditions, we expect a shift to forests dominated by species with relatively conservative traits, together with an associated decrease in forest annual increment, representing a challenge for future forest management.

| BEF relationships across different bioclimatic regions
We found substantial differences in the mechanisms by which functional diversity influences annual increment between the climatic regions. We found a more predominant effect of niche complementarity (i.e., functional diversity) on annual increment in the Mediterranean climate region. In other words, competitively dominant species (i.e., mass ratio effects) appear to be less important to forest functioning under harsh conditions (i.e., increased water limitation), which is consistent with the prediction of the stress gradient hypothesis in that the frequency of competitive interactions will vary inversely across abiotic stress gradients (Paquette & Messier, 2011;Rita & Borghetti, 2019;Wang et al., 2019).
The fact that we found a stronger effect of niche complementarity (i.e., functional diversity) on annual increment in water-limited plots could inform future forest management aiming to maintain annual increment under increasing drought (Spinoni et al., 2018). We showed that under increased water limitation, functional composition shifted to more conservative resource strategies, suggesting a decrease in forest annual increment. However, forest annual increment only showed a weak negative correlation with water limitation, shedding light on the importance of functional diversity for future forest management to maintain forest annual increment.

RM was supported by grant 2019-03758 from the Swedish Research
Council, Vetenskapsrådet.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The corresponding author confirms on behalf of all authors that there have been no involvements that might raise the question of bias in the work reported or in the conclusions, implications, or opinions stated.

FU N D I N G I N FO R M ATI O N
RM was supported by funding from the Swedish Research Council, Vetenskapsrådet (grant 2019-03758).

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in our manuscript are already hosted in publicly available archives. Specifically, the Italian National Forest Inventory data is available on the INFC website (https://www.inven tario fores tale. org/en). Functional trait data used in the study was downloaded from publicly available sources TRY (Kattge et al., 2020) and WOODIV (Monnet et al., 2021). Additionally, the r-code that support the findings of this study are publicly available at https://github.com/bobmu scare lla/INFC-funct ional -diver sity.